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1.  Introduction 


A  pacing  item  in  the  development  of  hypersonic  vehicle  design  is  the  capability  to  accurately 
predict  the  heating  on  cowl-lips  due  to  shock-on-shock  interactions  [1],  The  forebodies  of 
vehicles  proposed  for  high  speeds  typically  act  as  a  compression  system  which  facilitates 
optimum  mass  flow  through  the  inlet.  Depending  upon  the  flight  condition,  the  resultant 
system  of  forebody  shocks  may  interact  with  the  bow  shock  due  to  the  cowl-lip  of  the 
inlet.  The  resulting  interaction,  often  denoted  “interference  pattern,”  can  lead  to  severe 
amplification  of  pressure  and  heat  transfer  on  the  cowl-lip.  This  may  have  catastrophic 
consequences  for  the  vehicle  unless  accounted  for  in  the  design. 

Shock-on-shock  interactions  in  the  presence  of  a  cowl-lip  can  result  in  a  range  of  patterns 
depending  primarily  on  the  location  of  the  intersection  between  the  impinging  and  cowl 
shocks.  Edney  [2]  classified  these  distinct  interactions  into  six  different  categories:  Types  I 
to  VI.  The  highest  surface  loading  occurs  in  the  Type  IV  interaction  shown  schematically  in 
Figure  1.  This  pattern  is  obtained  when  the  impinging  shock  intersects  the  nearly  normal 
portion  of  the  cowl  bow  shock.  Two  triple  points  are  formed  from  which  emanate  two 
shear  layers.  The  flow  between  these  passes  through  a  series  of  weak  oblique  shocks  and 
expansion  fans  until  the  supersonic  jet  so  formed  strikes  the  surface  of  cylinder  (or  cowl-lip) 
after  the  formation  of  a  terminating  jet  bow  shock.  The  flow  then  expands  on  either  side  of 
the  stagnation  point.  The  jet  impingement  process  results  in  the  maximum  augmentation 
of  peak  heat  transfer  and  surface  pressure. 

Several  studies  have  investigated  shock-on-shock  interactions  at  a  range  of  free  stream 
Mach  numbers.  Following  the  efforts  of  Edney,  the  problem  has  been  addressed  with 
both  experimental  [■? - 5]  and  theoretical  efforts.  Among  the  latter  are  the  semi-empirical 
methods  of  Hains  et  al.  [(>],  Morris  et  al.  [7j  and  the  graphical  methods  of  Crawford  [8]  and 
B  ram  let  to  [9], 


Distorted  Bow  Shock 


Tx  =  200.8/? 
TwalI  -  530 R 
Impinging  shock  angle  =  18.1° 


Figure  1:  Schematic  of  Type  [V  interference  pattern 


Computational  solutions  with  the  Navier-Stokes  equations  have  been  reported  by  several 
researchers,  including  Tannehill  et  al.  [10],  White  et  al.  [11],  Klopfer  et  al.  [12],  Peery  et 
al.  [13],  Moon  et  al.  [11]  and  Thareja  et  al.  [15].  These  studies  employ  a  range  of  numerical 
algorithms  including  finite-difference  and  finite- volume  implementations  of  centered  and 
characteristic-based  schemes.  In  addition  to  providing  a  quantitative  assessment  of  the 
amplification  in  surface  loading,  these  solutions  have  advanced  the  understanding  of  the 
characteristics  of  shock-on-shock  interactions.  The  resulting  pattern,  for  example,  has  been 
shown  to  be  very  sensitive  to  the  location  of  the  point  of  intersection  of  the  impinging  and 
bow  shocks. 

One  of  the  few  unsteady  calculations  to  date  was  reported  by  Klopfer  et  al.  who 
examined  the  effect  of  varying  the  impinging  shock  location  dynamically,  as  would  be 
expected  in  a  flight  vehic  le.  As  the  impinging  shock  is  moved  from  one  extreme  position  to 
the  other,  all  of  the  patterns  of  Fxlney  were  successively  obtained.  Nevertheless,  a  common 
underlying  assumption  in  each  of  the  above  calculations  has  been  that  for  a  given  loc  ation 
of  the  impinging  shock,  i.c.,  for  a  given  interaction  pattern,  the  flow  is  steady. 

However,  there  is  reason  to  believe  that  some  Type  IV  interactions  are  unsteady.  A 
concrete  example  is  presented  in  Figure  1  which  focuses  on  the  Mach  8  Type  IV  interaction 
examined  in  this  research.  The  impinging  shock  is  created  upstream  by  a  wedge  of  15°  angle 
resulting  in  post-impinging  shock  flow  of  Mach  5.25.  The  cowl-lip  is  adequately  simulated 
by  a  cylinder.  This  configuration  was  investigated  experimentally  by  Wieting  et  al.  [5]  and 
has  been  examined  in  several  computational  studies  in  recent  years. 

Although  the  flow  was  assumed  to  be  steady  initially,  Holden  [16]  suggests  that  unsteadi¬ 
ness  in  the  frequency  range  between  I)  to  10kHz  is  evident,  in  surface  loading  experimental 
data.  Indeed,  several  computational  studies  have  noted  a  difficulty  in  obtaining  converged 
results  for  Type  IV  interac  tions  not  only  under  the  above  conditions  [12.13]  but  at  higher 
Mac  h  numbers  [17]  as  well.  However,  no  linkage  was  proposed  between  the  unsteadiness  in 


aggr^te  quantities  and  coherent  flow  structures. 

t  he  first  attempt  at  analyzing  the  details  of  the  unsteady  flow  was  performed  by  Kroll 
et  al.  [18].  Their  research  examined  the  interaction  under  inviscid  conditions  with  a  range 
of  upwind  and  centered  schemes.  In  a  detailed  study  that  included  several  time-accurate 
integration  methods  and  degrees  of  mesh  resolution,  they  simulated  the  unsteadiness  of  the 
interaction  and  related  it  to  the  coherent  flow  features.  They  concluded  that  the  unsteady 
limit  cycle  consists  of  an  oscillation  of  the  terminating  jet  how  shock  about  a  mean  posit  ion. 
This  is  accompanied  by  attendant  effects  in  the  remainder  of  the  flow.  Further  discussion 
of  their  results  is  presented  in  Chapter  4. 

This  inviscid  analysis,  while  useful  in  characterizing  the  unsteadiness,  fails  to  provide 
critically  important  information  about,  thermal  loading.  Further,  viscous  effects  aio  antici¬ 
pated  to  be  significant  in  the  stagnation  region  in  the  close  vicinity  of  the  terminating  jet 
bow-  shock.  The  impact  of  viscosity  on  the  nature  of  the  unsteadiness  is  not  immediately 
evident.  It  is  the  principal  purpose  of  this  research  to  address  these  questions  through 
numerical  simulation. 

In  recent  years,  significant  advances  iiave  been  attained  in  the  development  of  numerical 
schemes  for  the  Filler  equations.  Unlike  the  centered  schemes  of  Refs,  [lb  21],  the  newer 
schemes  use  upwind  met  hods  to  obviate  the  need  for  explicit  damping  terms.  These  newer 
methods  may  be  classified  into  the  means  by  which  higher  order  accuracy  is  attained 
MIJSCL  (Monotonic  Upwind- (’entered  Scheme  for  Conservation  Laws)  or  non-MUSCL. 
Examples  of  the  former  are  the  flux-vector  split  methods  of  Steger  and  Wanning  [22]  and 
of  van  Leer  [24]  and  the  flux-difference  split  method  of  Roe  [24],  Methods  obtaining  higher 
order  accuracy  with  non-MUSCL  methods  include  the  “upwind”  and  “symmetric"  TVD 
(Total  Variation  Diminishing)  methods  of  Yee  [25]  and  Marten  and  Yee  [26],  respectively. 

Although  preliminary  calculations  with  the  van  Leer  and  Roe  schemes  are  reported  in 
this  work,  the  principal  method  employed  is  that  of  MacCormack  and  Candler  [27],  denoted 


as  the  MC  method.  A  brief  description  is  provided  in  Chapter  2  followed  by  a  note  on 
the  boundary  conditions  and  other  numerical  details  (Chapter  d).  Finally,  the  results  are 
examined  in  Chapter  I  which  includes  an  investigation  of  the  effect  of  mesh  resolution  as 
well  as  of  the  time-accuracy  of  the  algorithm. 


2.  Theoretical  Model 


The  ‘2-D  full  Navier  Stokes  equations  are  solved  in  strong  conservation  form: 


dU  dF  d& 
Ot  di  +  drj 


=  0 


(2.1) 


The  £  coordinate  lines  wrap  around  the  cylinder  while  the  rj  lines  are  body-normal  thus 
forming  an  orthogonal  mesh.  In  Equation  2.1,  V  is  the  solution  vector  of  conserved  variables 
( U  =  { p ,  pu,  pv,  pt}j ./)  while  the  flux  vectors,  F  and  G  include  the  full  viscous  and  inviscid 
terms.  The  perfect  gas  assumption  is  invoked,  the  molecular  viscosity  is  approximated  by 
Sutherland’s  law  and  the  molecular  Frandtl  number  Pr  is  assumed  to  be  constant  (0.72). 

Equation  2.1  may  be  interpreted  as  describing  the  balance  of  mass,  momentum  and 
energy  fluxes  over  an  arbitrary  control  volume.  With  first-order  Euler  implicit  time  dis¬ 
cretization  and  linearization  of  the  fluxes  in  time,  the  discretized  equation  may  be  written 
in  the  form: 


{NUMERIC'S}  6U  =  PHYSICS  (2.2) 


where  PHYSICS  represents  the  residual,  NUMERICS  contains  the  driving  terms  and  bU 
denotes  the  change  in  the  solution  vector  at  each  time  step.  The  full  Navier-Stokes  equations 
are  utilized  in  computing  the  residual  with  the  appropriate  upwind  model  for  the  inviscid 
terms.  Viscous  flux  terms  are  centrally  evaluated. 

For  simplicity,  two  assumptions  are  employed  in  evaluating  the  NUMERICS  port  ion  of 
Equation  2.2.  First,  the  spatial  differences  in  flux  Jacobians  are  calculated  with  the  first 
order  spatially  accurate  Steger  Warming  method.  This  enhances  the  diagonal  dominance 
property  of  Equation  2.2  [28]  leading  to  a  robust  solution  procedure  as  described  below. 
Secondly,  the  thin  layer  approximation  is  employed  for  the  viscous  Jacobians  resulting,  in 
a  significant  reduction  of  computational  requirements. 

Time  integration  is  achieved  with  a  residual  driven  line  Causs-Seidel  relaxation  scheme 


as  described  in  Ref.  [29].  With  this  choice,  Equation  2.2  represents  a  block  tridiagonal 
system.  The  method  is  unconditionally  stable  for  mode!  diffusion  and  convection  equations 
and  is  also  known  to  be  relatively  insensitive  to  the  choice  of  time  increment  per  itera 
tion  [30].  The  formal  time-accuracy  (first  order)  of  the  resultant  algorithm  is  not  affected 
by  the  simplifications  described  above.  Nevertheless,  the  use  of  inconsistent  differencing 
needs  further  study  particularly  since  the  NUMERICS  portion  employs  first  order  spatial 
derivatives.  4  detailed  analysis  of  the  time  accuracy  of  the  method  is  deferred  to  Chapter  4. 

A  brief  description  of  the  spatial  scheme  is  presented  with  reference  to  the  evaluation 
of  the  ('artesian  flux  (G)at  a  j  -f  ^  surface  of  the  cell  centered  finite  volume  formulation. 
The  extension  to  general  coordinates  is  straightforward  [31]  as  is  the  implementation  of 
the  viscous  terms  which  are  centered.  At  each  interface,  the  state  of  the  flow  is  described 
by  two  vectors  of  conserved  variables,  UL  and  U on  either  side  of  the  interface.  These 
are  obtained  by  upwind  extrapolation  to  the  cell  surface  with  the  MUSCL  approach  of 
van  Leer  [32]  in  conjunction  with  a  limiter.  With  W  denoting  the  vector  of  extrapolated 
quantities: 


Ki  = 

where  A  introduces  the  limiter  function,  /: 


"  2AW 


j  —  /(Aj+l,  A,_ i);  AJ+i  -  WJ+i  Wj ;+l 


i  '  3) 


(2,  , 


In  the  present  work  the  primitive  variables  (W  =  {p, u,v,p})  are  extrapolated  and  U 
is  then  reconstructed.  The  limiter  /,  preserves  monotonicity  within  the  solution.  Many 
limiter  functions  have  been  proposed  in  the  literature.  The  calculations  in  this  work  employ 
the  minmod  limiter  [33]  chosen  for  its  popularity  and  robustness: 


f{x,y)  -  minmod  (or,  y)  =  sgn  (x)  •  max  {0,min[|j|,j/  sgn(.r)]}  (2.5) 


Following  Steger  et  al.  [22],  the  first  order  homogeneous  hyperbolic  property  of  the 
flux  Jacobian  B ,  ( B  =  dG/dU,  G  =  BU ),  is  utilized  to  split,  the  fluxes  into  positive  and 


i 


negative  components: 


B  =  Q~lAQ  =  Q-‘{  A+  +  A_)Q  =  D+  +  i?_  (2.6) 

where  A  is  a  diagonal  matrix  consisting  of  the  eigenvalues  of  B,  A  =  diag  {i>,r  +  c,  e  -  r}. 
A+  and  A_  denote  the  splitting  of  the  eigenvalues  into  positive  and  negative  components. 
At  a  face  j  +  1/2,  therefore,  the  original  Steger  Warming  (abbreviated  SW)  flux  {Gsw) 
may  be  written  as: 

6V,+l  =  Bl;VL  +  B*UR  (2.7) 

The  inviscid  flux  formula  of  MacCormack  and  Candler  [27]  (Gmc) 

LiM.  ,  tifl  „ 

Gmoj+j  =BS  U'-  +  B.’  V"  (2.8) 

where  UL  and  UR  are  evaluated  with  the  MUSCL  approximation  (Eqns  2.3  and  2.4).  The 
essential  difference  between  the  SW  and  the  MC  schemes  is  in  the  manner  in  which  the 
flux  Jacobians  ( B )  are  evaluated.  While  SW  employs  upwinding  for  the  Jacobians,  MC  is 
centered.  Additional  details  of  the  MC  method  may  be  found  in  Ref.  [34]. 

The  MC  scheme  was  developed  for  application  only  in  the  boundary  layers,  and  in  fact 
leads  to  instability  when  applied  in  regions  of  discontinuities  such  as  shock  waves.  As  a 
result,  it  is  necessary  to  revert  to  the  SW  scheme  near  shocks.  This  is  achieved  in  a  smooth 
manner  by  defining  a  pressure  based  switch, 

Xt  =  1  +  PRxPR  (2-9) 

PR  =  -  (2.10) 

and  defining  the  flux  at  the  j  +  1/2  face  as: 

Gj+i  =  Xi  Gmc  +  (1  —  \i  )G\sn'  (2.11) 
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3.  Boundary  Conditions  and  Numerical  Details 


The  boundaries  may  be  categorized  into  the  following.  At  an  inflow  boundary,  the  flow 
vector  {p,  pu,pv,pe]  is  specified  according  to  the  known  freestream  values.  At  solid  bound¬ 
aries,  the  velocity  vector  and  the  normal  pressure  gradient  are  assumed  to  be  zero  and  at  a 
constant  surface  temperature.  At  outflow  boundaries,  the  flow  is  assumed  predominantly 
supersonic  and  the  zero  gradient  extrapolation  condition  is  applied.  The  boundary  condi¬ 
tions  for  the  implicit  portion  of  the  algorithm  are  described  in  Ref.  [35]. 

The  convergence  or  the  absence  of  it  for  unsteady  calculations  is  determined  by  moni¬ 
toring  several  quantities.  The  global  residual,  defined  as: 


\\G.R.\\  = 


1 


IL  JL  4 


i  jt=i  \  Un.k 


(3.1) 


is  a  measure  representative  of  the  overall  solution.  In  Equation  3.1,  k  denotes  the  kth 
component  of  U  on  an  IL  x  JL  computational  mesh  and  (Rk)i,j  is  the  change  in  the  flow 
solution,  Rk  =  at  i,j  and  —  {p,/tw,pu,pe ■}0O.  For  plotting  purposes,  these  values 

are  normalized  by  the  value  of  j|6'.i2.||  obtained  after  the  first  iteration.  Convergence  is 
assumed  when  ||6\i?.||  drops  S  or  more  orders  of  magnitude.  After  this,  the  surface  pressure 
and  heat  transfer  are  monitored  over  several  two  characteristic  times,  Tc  =  L / U<x, ,  where 
Uoo  is  the  freestream  velocity  and  L  is  a  macroscopic  length  scale  -  the  diameter  of  the 
cylinder  or  the  length  from  the  leading  edge  to  the  corner.  In  addition,  the  integrated 
root  mean  square  (RMS)  pressure  and  heat  transfer  values  over  the  entire  surface  are  also 
monitored.  With  P  and  Q  denoting  pressure  and  heat  transfer  respectively: 
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4.  Results 


The  parameters  for  this  flow  are  chosen  after  consideration  of  available  experimental  data 
from  Wieting  et  al.  [5]:  =  8.03,  Cylinder  Radius  =  1.5in.,  =  200.8/?  and  Twa||  - 

530 R.  The  impinging  shock  angle  is  18.1°  corresponding  to  a  post-shock  Mach  number  of 
5.25.  The  experimental  results  were  reanalyzed  and  published  by  Thareja  e.t  al.  [15]  who 
also  performed  a  computational  study  in  which  they  varied  the  impinging  shock  location 
to  obtain  the  best  fit  with  experimental  data.  From  their  work,  we  choose  the  reference 
point  for  the  impinging  shock  as  (— 3.5/ 1 .5/2,  -0.522/1.5/2)  where  R  is  the  radius  and  the 
point  (0,0)  corresponds  to  the  cylinder  center.  This  is  denoted  as  Case  SO  in  Ref.  [15]. 

Figure  2  illustrates  the  general  features  of  the  mesh  employed.  The  inflow  computa¬ 
tional  boundary  is  assumed  to  be  a  sequence  of  piecewise  smooth  spirals  bounding  the 
domain  of  interest.  In  the  radial  direction  the  grids  are  stretched  algebraically.  Three  grids 
are  employed  in  this  study  with  the  coarser  meshes  (Grids  1  and  2)  being  extracted  sys¬ 
tematically  from  the  finest  mesh  (Grid  3)  in  a  multigrid  fashion.  The  salient  characteristics 
of  each  grid  are  presented  in  Table  1. 

In  the  coarse  mesh  calculation  the  basic  features  of  the  flow  structure  are  reproduced 
(Figure  3).  However,  the  shock  system  exhibits  significant  staircasing  and  the  details 
of  the  supersonic  jet  and  the  bounding  shear  layers  are  not  discernible.  One  interesting 
aspect  of  the  calculation  is  the  dependence  of  the  final  solution  upon  the  time  step  size 
utilized.  At  low  time  steps,  an  unsteadiness  is  evident  in  the  global  residual  criterion  and, 
to  a  lesser  extent,  in  the  aggregate  surface  quantities  (RMS  pressure  and  heat  transfer). 
On  the  other  hand,  at  larger  time  step  sizes,  the  solution  converges  to  machine  accuracy. 
Preliminary  calculations  with  both  the  van  Leer  and  Roe  methods  exhibit  unsteady  limit 
cycles  in  each  of  the  monitored  quantities  at  all  achievable  time-step  sizes.  The  employment 
of  characteristic  boundary  conditions  on  the  outflow  boundaries  has  no  effect  upon  these 
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Table  1:  Shock  Interference  Grid  Details 


Grid 

IL  x  JL 

R^c\max 

|au 

a  e\m,n 

A0]mai. 

EB 

1 

49  x  24 

5.6 

17.9 

10.4 

1.8 

4.0 

3.1 

2.0 

2 

99  x  49 

2.2 

KiH 

4.3 

0.9 

2.0 

1.5 

1.4 

3 

199  x  99 

1.0 

3.4 

2.0 

0.45 

1.0 

0.8 

1.2 

Legend:  IL  -  Points  in  £  direction  JL  -  Points  in  ij  direction 

Rec  -  Surface  mesh  Reynolds  number  A 8  -  Angular  spacing  (deg) 
c  -  Stretch  factor  at  surface 

Subscripts:  av  -  average  min  -  minimum 

max  -  maximum 

limit  cycles  which  arc  described  in  more  detail  below. 

The  apparent  contradiction  in  the  asymptotic  behavior  among  the  three  algorithms 
(MG.  van  Leer  and  Roe)  is  examined  further  with  reference  to  the  numerical  method. 
For  problems  with  steady-state  solutions,  the  code  is  typically  operated  at  the  highest 
achievable  inviscid  OFL  number.  The  MC  scheme  is  most  robust  in  this  regard.  Operational 
CFL  numbers  between  1000  —  5000  are  achieved  with  the  higher  number  pertaining  to  the 
coarser  meshes.  In  contrast,  achievable  CFL  numbers  with  Roe’s  scheme  are  an  order  of 
magnitude  lower  while  the  van  Leer  scheme  is  the  least  robust  ( CFL  <  25).  It  is  known 
that  at  high  time-step  sizes  and  with  a  fixed  number  of  sweeps,  the  Gauss-Seidel  algorithm 
exhibits  maximal  damping  [30]  indicating  the  possibility  of  false  steady  results  with  higher 
A t  values.  On  the  other  hand,  Yee  et  al.  [36]  state  that  “discrete  maps”  of  nonlinear 
equations  can  produce  solutions  which  do  not  satisfy  the  original  continuous  system.  Such 
solutions  can  manifest  themselves  as  spurious,  unsteady  cycles  of  any  period,  and  may  be 


obtained  at  time  step  sizes  above  or  sometimes  even  below  the  linearized  stability  limit. 
They  suggest  that  unsteady  asymptotic  behavior  be  examined  in  the  limit  n  — ♦  oo,  A#  — +  0 
where  n  is  the  number  of  interations. 

The  time  accuracy  of  the  present  code  is  first  investigated  on  the  intermediate  mesh.  The 
discretization  utilized  is  based  on  the  Euler  implicit  formula  and  is,  therefore,  formally  first- 
order  time  accurate.  The  solution  procedure,  CJauss-Seidel  line  relaxation,  may  be  viewed 
as  time  accurate  within  the  framework  of  the  sub-iteration  strategy  [28].  Nevertheless.  I  he 
use  of  inconsistent  differencing  (through  the  use  of  Sieger  Warming  Jacobiaus)  on  the  left 
hand  NUMERICS  portion  is  relevant  to  this  issue  anti  a  detailed  examination  is  performed 
in  the  limit  At  — >  0  by  a  combination  of  increased  number  of  sweeps  per  iteration  (from  two 
to  four)  and  a  reduction  in  At.  At  each  time-step,  therefore,  the  equations  are  solved  to 
within  a  certain  tolerance.  A  consistent  strategy  of  reducing  At  is  followed  until  the  results, 
monitored  with  the  above  mentioned  aggregate  and  local  quantities,  are  independent  of  At. 
The  assumption  that  the  limit  cycle  is  independent  of  the  initial  conditions  is  implied  and 
tested  by  starting  the  soli,  'on  from  a  converged  blunt  body  shock  and  separately  from  a 
uniform  flow  condition.  The  time-step  size  independent  limit  cycles  in  both  instances  are 
identical. 

The  process  outlined  above  is  discussed  with  reference  to  Grid  2.  A  limit  cycle  is  first 
obtained  with  a  relatively  high  but  constant  At,  corresponding  to  a  CFL  number  varying 
slightly  about  4000.  The  time-step  is  successively  halved  and  the  solution  marched  until 
a  fiew  limit  cycle  is  reached.  This  process  is  halted  when  halving  At  has  no  perceptible 
influence  upon  the  limit  cycle  as  reflected  on  several  aggregate  and  local  quantities.  The 
results  are  shown  in  Figure  4(a)  for  the  quantity  log(||Gil||)  of  Equation  3.1.  The  extreme 
values  in  the  residual  exhibit  a  four-period  type  of  oscillation.  Although  the  profiles  of  the 
other  quantities  monitored  are  different,  the  influence  of  time-step  reduction  is  similar.  The 
initial  reductions  in  At  result  in  an  increase  in  amplitude  of  oscillation.  Below  At  =  0.0032).. 


14 


corresponding  to  a  CFL  number  of  about  125,  no  further  influence  of  reduction  is  observed. 
The  implicit  solution  is  compared  with  results  from  an  explicit  second  order  time  and  span* 
accurate  method  -  a  two-step  predictor  corrector  algorithm  (Heim’s  scheme).  This  method 
yields  time-step  size  independent  results  at  the  linearized  stability  limit  (CFL  ~  0.9).  The 
entire  log(||G.R||)  history  of  the  explicit  algorithm  is  shown  in  Figure  4(b). 

Results  from  the  explicit  and  implicit  algorithm  are  compared  in  Figure  5  for  several 
quantities  monitored  over  one  cycle  spanning  the  dominant  frequency.  In  these  figures,  the 
profiles  from  the  explicit  calculation  are  translated  along  the  time  axis  until  they  match  the 
implicit  time  development.  The  displayed  quantities  are  the  aggregate  quantities:  global 
residual,  the  instantaneous  RMS  surface  pressure  and  heat  t  ransfer  over  the  entire  surface 
and  the  local  quantities:  maximum  surface  pressure  and  heat  transfer  which  occur  in 
the  vicinity  of  the  location  of  the  impinging  jet.  It  is  clear  that  the  error  in  amplitude 
and  frequency  between  the  implicit  and  explicit  methods  is  negligible  for  all  quantities 
monitored.  Indeed,  the  largest  discrepancy  occurs  in  the  maximum  heat  transfer  rate  on  t  he 
surface  where  the  relative  error  is  only  0.1%.  Although  the  unsteadiness  in  this  computation 
is  evident  in  the  aggregate  quantities  presented,  the  variation  of  local  quantities  is  not 
pronounced.  For  example,  the  difference  between  the  peak  and  trough  in  the  maximum 
heat  transfer  cycle  is  only  2%  of  the  mean  while  that  for  pressure  is  1.25%.  The  detailed 
flow  features  reveal  only  a  modest  degree  of  unsteadiness  in  the  overall  flow  structure. 

The  unsteady  features  are  more  apparent  in  results  obtained  on  Grid  3.  These  re¬ 
sults  were  obtained  at  At  =  0.00lTf.  Figure  6  displays  the  computed  limit  cycles  of  the 
log(||6'jR||)  and  maximum  pressure  and  maximum  heat  transfer,  respectively.  The  basic 
four-period  type  oscillation  similar  to  that  observed  on  Grid  2  is  discernible  although  the 
regularity  does  not  persist.  The  range  of  values  obtained  by  the  maximum  heat  transfer 
and  pressure  is  much  larger,  however  (65%  and  30%  about  the  mean,  respectively).  A 
spectral  analysis  of  these  limit  cycles  reveals  a  dominant  frequency  of  about  32kHz.  In 
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comparison,  the  experimentally  observed  frequencies  for  Typo  IV  interactions  are  in  the 
range  of  3  to  tOkllz  [16].  An  estimate  of  the  frequency  supportable  by  Mesh  3  may  be 
based  on  a  wavelength  of  2  Ax  and  a  sonic  wave  speed  in  the  stagnation  region,  (’boosing 
A.r  as  the  grid  spacing  in  the  circumferential  direction,  the  resolvable  frequency  ( 1080k Hz ) 
is  much  higher  than  that  observed  in  the  present  physical  phenomenon.  The  number  of 
time-steps  within  each  cycle  is  about  800. 

A  small  portion  of  the  residual  cycle  spanning  two  consecutive  crests  of  maximum 
pressure  representing  the  dominant  frequency  is  isolated  for  further  study  (Figure  7(a)). 
The  maximum  heat  transfer  cycle  is  shown  in  Figure  7(b).  The  time  between  the  two 
pressure  peaks  is  0.75Tr  corresponding  to  a  frequency  of  The  maximum  heat 

transfer  is  similar  to  the  peak  pressure  with  a  small  leading  phase. 

Five  points  denoted  A,B,C,D  and  E  are  identified  on  the  pressure  profile  between  the 
two  successive  peaks.  The  Mach  contours  at  these  points  are  shown  in  Figure  8.  The 
distorted  bow  shock  as  well  as  the  impinging  shock  are  captured  within  at  most  two  zones. 
The  shear  layers  emanating  from  the  two  triple  points  are  also  well  resolved.  However,  some 
portions  of  the  distorted  bow  shock  exhibit  the  “zig-zag"  phenomenon  which  can  result  in 
a  degradation  of  accuracy  [37]. 

The  most  significant,  difference  in  the  Mach  patterns  appears  to  be  the  angle  made  by 
the  terminating  jet  shock  with  the  surface  of  the  cylinder.  This  shock  oscillates  about  two 
positions.  At  one  extreme  (location  zl)  it  forms  an  acute  angle  with  the  surface  tangent  of 
I  he  cylinder  at  the  stagnation  point.  At  the  other  extreme  (location  C),  the  terminating  jet 
shock  is  parallel  to  the  cylinder  surface  tangent.  This  shock  movement  is  accompanied  by 
changes  in  the  angle  with  which  the  supersonic  jet  between  the  two  shear  layers  impinges 
upon  the  body  surface. 

The  surface  pressure  and  beat  transfer  corresponding  to  the  points  .4  through  E  are 
exhibited  in  Figure!).  The  picture  is  normalized  by  the  stagnation  pressure  obtained  with 
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Figure  M:  Mad)  contours  at  various  points  in  limit  cycle  •  Grid  il  Type*  IV 
interaction 
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Figure  9:  Surface  loading  at  various  points  in  limit  cycle  -  Grid  l 
interaction 


the  Uayleigh  formula.  For  heat  transfer  normalization,  we  follow  the  approach  of  Thareja 
fit  al.  [15].  Although  the  experimental  noninterfering  stagnation  heat  transfer  value  is 
61.7 Dtu/ft2  —  sec,  they  argue  that  since  the  theoretical  value  for  peak  heat  transfer, 
obtained  with  a  VSL  (Viscous  Shear  Layer)  analysis  is  MAWBtuj ft1  —  .sec,  the  difference 
may  be  factored  out  by  utilizing  the  higher  (61.7)  value  to  normalize  experimental  results 
and  the  lower  value  (41.43)  for  computed  values.  The  abscissa,  0.  is  the  angle  measured 
with  respect  to  the  horizontal,  negative  values  denoting  the  cylinder  surface  below  the 
noninterfering  stagnation  streamline  (Figure  8).  At  point  .4,  the  first  peak  of  the  pressure 
cycle,  the  terminating  jet  points  upward  and  the  corresponding  maximum  pressure  occurs  at 
9  ~  —18.5°.  Between  points  A  and  B ,  the  terminating  jet  shock  assumes  a  posit  ion  parallel 
to  the  surface  tangent  and  remains  so  until  after  point  C.  During  this  time,  the  maximum 
surface  pressure  location  moves  downward  to  a  position  of  6  ~  —21°  and  the  terminating 
jet  impinges  more  perpendicularly  to  the  surface  than  at  point  .4.  Between  points  C  and  D. 
the  terminating  jet  shock  again  assumes  a  position  similar  to  that  at  point  .4  and  remains 
in  this  position  until  point  E.  Correspondingly  the  maximum  surface  pressure  point  moves 
upwards  until  it  reaches  a  maximum  of  —18.2°.  The  temporal  development  of  the  surface 
heat  transfer  is  generally  similar  to  that  of  the  surface  pressure  in  phase  and  location.  After 
point  £\  the  cycle  repeats  in  t  he  qualitative  features  though  the  peak  values  obtained  differ 
by  a  modest  amount  between  various  cycles. 

Overall,  the  computed  surface  pressure  over  predicts  the  experimental  observations.  The 
computed  location  of  the  maximum  pressure  oscillates  A 6  ~  2.5”  about  a  mean  at  9  ~ 
—20°.  The  best  agreement  in  magnitude  of  the  pressure  peak  is  obtained  at  point  C  on 
the  limit  cycle  (6%  over  predicted).  However,  the  location  of  this  peak  (  —  19.5°)  is  about 
4.5”  above  the  experimental  value  (—24'').  The  maximum  discrepancy  in  peak  pressure  is 
at  point  E. 

In  contrast  to  the  surface  pressure,  the  heat,  transfer  is  underpredicted  at  all  points  in 


the  limit  cycle.  One  explanation  is  based  on  the  fact  that  the  impinging  supersonic  jet  is 
processed  through  several  oblique  shocks  and  expansion  fans.  In  this  region,  the  effect  of 
the  switching  mechanism  reverting  the  MC  scheme  to  the  SW  scheme  is  not  clear  and  may 
influence  the  quantitative  calculations  in  the  surface  quantities.  Inadequate  grid  resolut  ion 
may  be  another  factor.  The  lowest  value  of  peak  heat  transfer  occurs  at  point  C'  while 
the  highest  value  (underpredicted  by  20%)  occurs  at  point  E.  The  locations  of  the  heat 
transfer  peaks  fall  roughly  within  the  experimental  envelope.  In  the  lower  portion  of  the 
cylinder,  in  the  region  between  —25°  and  —45°,  both  the  surface  pressure  and  heat  transfer 
distributions  exhibit  spatial  as  well  as  temporal  oscillations.  The  heat  transfer  in  this  region 
as  well  as  on  the  upper  portion  of  the  cylinder  is  underpredicted.  This  is  similar  to  the 
overall  observations  of  other  researchers  [12]. 

Further  examination  of  the  computed  flowfield  reveals  the  existence  of  two  distinct  re¬ 
gions  of  unsteady  separation,  one  below  the  jet  impingement  position  (—35°  to  -  29") 
and  the  other  above  (24°  to  35°).  The  lower  separated  flow  is  in  the  vicinity  of  the  region 
where  the  surface  pressure  and  heat  transfer  were  seen  to  exhibit  spatial  oscillations.  The 
angular  extent  of  this  separation  region,  monitored  with  the  skin  friction  coefficient,  fluc¬ 
tuates  in  time.  The  second  separation  region,  above  the  stagnation  point,  is  caused  by  a 
recompression  experienced  by  the  flow  as  it  rounds  the  upper  portion  of  the  cylinder.  This 
compression  results  in  a  pressure  rise  on  the  surface  causing  boundary  layer  separation. 

It  should  be  noted  that  in  previous  research  with  the  Kuler  equations  [18],  a  recompres¬ 
sion  shock  system  was  observed  in  clear  detail  on  the  upper  portion  of  the  cylinder.  'The 
surface  Mach  number  near  the  point  of  impingement  of  the  compression  shock  was  higher 
than  unity.  This  is  precluded  in  the  present  research  by  the  viscous  boundary  conditions. 
The  inviscid  results  of  Ref.  [18]  for  the  same  interaction  with  the  Roe  as  well  as  with  the 
“Upwind”  TVD  schemes  also  indicate  an  oscillating  surface  pressure  profile.  In  that  calcu¬ 
lation,  however,  the  point  in  the  cycle  where  the  jet  terminates  normal  to  the  surface  (e.g., 
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point  C  above)  results  in  the  highest  pressure  peak  in  contrast  to  the  present  behavior. 
One  possible  explanation  for  this  difference  is  based  upon  the  influence  of  dissipation  upon 
the  strength  and  location  of  the  terminating  jet  shock.  At  point  C  in  the  limit  cycle,  the 
jet  shock  is  nearly  parallel  to  and  farther  away  from  the  surface  than  at  point  A  where  the 
shock  is  somewhat  oblique  to  the  flow.  In  the  viscous  stagnation  region,  the  dissipation 
mechanism  (whether  numerical  or  physical)  not  only  alters  the  total  pressure  and  tempera¬ 
ture  recovery  based  upon  shock  strength  but  may  also  change  the  information  propagation 
characteristics.  This  mechanism  could  introduce  a  phase  lag  between  the  shock  system  and 
the  corresponding  surface  effects. 

The  source  of  the  unsteadiness  observed  in  the  computation  is  not  clear.  The  numerical 
aspect  of  the  solution  has  not  been  fully  examined  in  the  context  of  spatial  and  temporal 
accuracy.  The  physical  mechanisms  acting  upon  the  triple-point  and  the  flow  properties  of 
the  shear  layers  bounding  the  supersonic  jet  are  complex.  In  addition  the  shear  layers  are 
subject  to  inviscid  instability  and,  depending  upon  the  Reynolds  number  in  the  adjacent 
shock  layer,  transition  to  turbulence  can  occur  [16].  Further  work  is  required  to  resolve 
these  issues. 


5.  Conclusions 


Results  for  the  Type  IV  interaction  display  unsteady  (low  behavior  with  all  three  algo¬ 
rithms.  A  study  of  the  computed  unsteady  flow  field  has  been  performed  with  the  MC 
scheme.  A  large  scale  movement  of  the  supersonic  impinging  jet  is  observed  with  accom¬ 
panying  influences  on  the  remainder  of  the  flow.  The  dominant  frequency  of  the*  observed 
phenomenon  is  about  32  kHz  on  the  finest  mesh.  Unsteady  separation  regions  are  ob¬ 
served  both  above  and  below  the  point  of  jet  impingement.  Surface  pressures  show  modes! 
discrepancies  with  experiment  while  heat  transfer  rates  are  underpredicted. 
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Nomenclature 

D  flux  Jacobian  of  G 

Dtu  British  Thermal  Unit 

c  stretch  factor;  local  speed  ol  sound 

cu  specific  heat  at  constant  volume 

C'FL  Courant- Friedrich- Levy  number 

c  total  energy 

( ,  internal  energy 

/  feet 

F, F,G,G  flux  vectors 

G. R.  global  residual 

IL,JL  points  in  £  and  //  direction 

J  Jacobian,  inverse  cell  volume 

A/  Mach  number 

MC  MaeCortnack  and  Candler  scheme 

SW  Original  Sieger  Warming  scheme 

p  pressure 

psi  pounds  per  square  inch 

Pr  Prandtl  number 

PR  pressure  function 

Q  matrix  diagonalizing  -4;  heat  transfer 

R  gas  constant;  Hankiue 

R<  Reynolds  number 

RMSP  root  mean  square  surface  pressure 


ill 


RMSQ 

root  mean  square  surface  heat  transfer 

t 

time 

T 

temperature 

Tc 

characteristic  time 

u 

Cartesian  velocity  in  x  direction 

u,u 

solution  vector 

V 

Cartesian  velocity  in  y  direction 

V 

Cartesian  velocity  vector 

VSL 

Viscous  Shear  Layer  analysis 

W 

vector  of  extrapolated  variables 

x,y 

Cartesian  coordinates 

6 

cutoff  value;  any  change 

A 

change  across  interface 

6 

cutoff  parameter 

V 

gradient  operator 

7 

ratio  of  specific  heats 

A 

eigenvalue 

A 

eigenvalue  matrix 

d 

partial  derivative  operator 

P 

density 

e 

angle 

transformed  coordinates 
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pressure  switch 


Xi 


Subscripts 

J  +  5 


w 


oo 


interface  between  cells  j  and  j  +  1 

viscous 

wall 

frcestream  conditions 


Superscripts 


L 

R 


state  at  left  of  interface 
state  at  right  of  interface 
positive  and  negative  components 
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